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^ Abstract 

^^ From a mathematical perspective, radiation hydrodynamics can be thought of as a 

CIh system of hyperbolic balance laws with dual multiscale behavior (multiscale behavior 

associated with the hyperbolic wave speeds as well as multiscale behavior associated 

with source term relaxation). With this outlook in mind, this paper presents a hybrid 



> 



^O Godunov method for one- dimensional radiation hydrodynamics that is uniformly well 

^Sl behaved from the photon free streaming (hyperbolic) limit through the weak equilib- 



in 

o 



rium diffusion (parabolic) limit and to the strong equilibrium diffusion (hyperbolic) 
limit. Moreover, one finds that the technique preserves certain asymptotic limits. The 
^^ method incorporates a backward Euler upwinding scheme for the radiation energy den- 

sity Er and flux F^ as well as a modified Godunov scheme for the material density />, 
momentum density m, and energy density E. 



'i 



The backward Euler upwinding scheme is first-order accurate and uses an implicit 
HLLE flux function to temporally advance the radiation components according to 
the material flow scale. The modified Godunov scheme is second-order accurate and 
directly couples stiff source term effects to the hyperbolic structure of the system of 
balance laws. This Godunov technique is composed of a predictor step that is based 
on Duhamel's principle and a corrector step that is based on Picard iteration. The 
Godunov scheme is explicit on the material flow scale but is unsplit and fully couples 
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matter and radiation without invoking a diffusion-type approximation for radiation 
hydrodynamics. This technique derives from earUer work by Miniati & Colella 2007. 
Numerical tests demonstrate that the method is stable, robust, and accurate across 
various parameter regimes. 



1 Introduction 

Radiation hydrodynamics is a dynamical description of fluid material interacting with elec- 
tromagnetic radiation and is appropriate whenever radiation governs the transport of energy 
and momentum in the fluid. Many phenomena in plasma physics and astrophysics are 
governed by radiation hydrodynamics, some examples include: star formation, supernovae, 
accretion disks, radiatively driven outflows, stellar convection, and inertial confinement fu- 
sion [5,37,40 . In these applications, the radiation field heterogeneously couples to the ma- 



terial dynamics such that radiative effects are strong in some parts of the system and weak 
in other parts. These variations give rise to characteristically different dynamical properties 
(i.e., advection versus diffusion behavior). The primary objective for developing the numer- 
ical technique presented in this paper is to have a computational tool that accurately solves 
radiation hydrodynamical problems across a range of asymptotic limits. The new algorith- 
mic ideas are cast in such a way that they seem familiar with respect to classical Godunov 
schemes and can be implemented in existing codes with minimal computational overhead. 
A future research endeavor is to combine the hybrid Godunov method for radiation hydro- 
dynamics with an existing code for MHD (magnetohydrodjTiamics) such as Athena 52 and 
investigate full radiation MHD in multiple spatial dimensions. 

Some of the initial research in developing numerical methods to solve radiation hydrodynami- 
cal problems was carried out by Castor 1972, Pomraning 1973, Levermore & Pomraning 1981, 
Mihalas Sz Klein 1982, and Mihalas & Mihalas 1984. One of the simplest and most successful 
approaches used in astrophysics was the Zeus code of Stone et al 1992, which relied on opera- 
tor splitting and a Crank- Nicholson temporal finite difference scheme. Since the introduction 
of that code, finite volume schemes (e.g., Godunov- type methods) have emerged as a pow- 
erful technique for solving hyperbolic conservation laws (i.e., the mathematical framework 



describing radiation hydrodynamics) [27,28 . Moreover, this integral formulation allows one 
to more naturally treat boundary conditions, capture shock waves and other discontinuous 
behavior, investigate complex geometries and multi-dimensions, and implement adaptive 
mesh refinement [59]. Despite these advantages, there have been significant difficulties in 
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developing a Godunov method to accurately represent radiation hydrodynamical behavior 
across a range of asymptotic limits. Earlier attempts to construct a Godunov-type technique 
either (i) neglected the heterogeneity of the matter-radiation coupling and solved the system 
of equations in a specific limit (9|[l0], (ii) were based on a model system [3||23], (iii) invoked 
artificial coupling terms that were based on the reference frame in which the problem was 
solved [1] as pointed out by Lowrie & Morel 2001, or (iv) used a variation of flux limited 
diffusion [22|[25|[29|[58] . Lowrie & Morel 2001 was even critical of the likelihood of developing 
such a method for full radiation hydrodynamics. 

Developing a Godunov method for radiation hydrodynamics has been difficult because nu- 
merical difficulties arise from (i) the reference frame one chooses for taking moments of the 
photon transport equation, (ii) multiscale waves, {in) stiff source terms, and (iv) solving a 
hierarchy of radiation transport moment equations to compute the variable tensor Edding- 
ton factor f . The first difficulty occurs because one takes moments of the photon transport 
equation in order to define the radiation quantities that interact with the material compo- 
nents of the system. One encounters problems with the radiation field's specific intensity 
function Ilu, n) which is governed by frame-dependent quantities, the radiation frequency u 
and directional vector n. Here, one either casts the photon transport equation into the co- 
moving frame (at rest with respect to the local fiuid velocity) and contends with complicated 
transport operators but simple interaction terms S{i>,n); or one casts the photon transport 
equation into the Eulerian frame (at rest with respect to the system as a whole) and contends 
with simple transport operators but complicated interaction terms. The specific intensity is 
written below in the Eulerian frame, where c is the speed of light and t is time: 

~ + n-v]liu,n) = Siu,n). (1) 

The Eulerian frame approach is also referred to as the mixed frame approach if the radiation 
intensity is measured in the Eulerian frame while the opacities a (embedded in the interac- 
tion terms) are measured in the comoving frame [6,31 . The second difficulty arises because 



although the transport operators have a simple formulation in the mixed/Eulerian frame, 
the radiation quantities are still characterized by waves propagating at the speed of light c. 
This dynamical scale is significantly larger than the speed of sound a^o which characterizes 
the material quantities in the absence of radiation, such variation in propagation speeds 
defines the nature of multiscale waves associated with radiation hydrodynamics. It is impor- 
tant to note that the reference frame one chooses to take moments of the photon transport 
equation does not affect how one defines the material quantities. A mixed frame approach 
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was adopted because the resulting equations most closely resemble a system of hyperbolic 
balance laws which is advantageous for constructing a Godunov-type method. The third dif- 
ficulty occurs because in addition to defining balance laws for the radiation quantities, one 
must also add relativistic stiff source terms that are correct to 0{u/c) to the right hand sides 
of the non-relativistic conservation laws for the material quantities (i.e., Euler equations). 
These source terms are stiff because of the variation in time and length scales associated with 



radiation hydrodynamical problems 38 . Having to contend with stiffness arising from some 
waves propagating at the speed of light as well as source terms having large magnitudes, it 
is obvious to see how such numerical difficulties make conventional techniques like operator 
splitting and the method of lines breakdown (27|[28|[56] . 

This paper presents a hybrid Godunov method that addresses the above numerical diffi- 
culties. The technique adopts a mixed frame approach, includes the appropriate frame- 
dependent terms to 0{u/c), is implicit with respect to the fastest hyperbolic wave speeds, 
and semi-implicitly updates the stiff source terms. The paper proceeds in the following 
manner. After defining the full system of equations for radiation hydrodynamics, the paper 
discusses what dynamics characterize the various asymptotic limits. Then, the paper gives 
an overview of the hybrid Godunov method and explains certain numerical properties that 
the algorithm possesses. The next three sections present the main algorithmic ideas behind 
the hybrid Godunov method. Lastly, the paper describes numerical tests that demonstrate 
the technique to be stable, robust, and accurate across various parameter regimes. 

2 Radiation Hydrodynamics 

As presented in Lowrie et al 1999 and Lowrie & Morel 2001, the system of equations for 
radiation hydrodynamics can be non-dimensionalized with respect to the material flow scale 
so that one can compare hydrodynamical and radiation effects as well as identify terms 



that are 0{u/c) [31,33 . This scaling gives two important parameters: C = c/a^ which 
measures relativistic effects and P = arT^/pooO^ which measures how radiation affects 
material dynamics. Additionally, a^ = Sir^k^/lbc^h^ = 7.57 x 10~^^ erg cm~^ K'"^ is a 
radiation constant. Too is a reference material temperature in the absence of radiation, and 
Poo is a reference material density in the absence of radiation. The full system of equations 
for radiation hydrodynamics in the mixed frame that is correct to 0{1/C) is: 

|^ + V-(m) = 0, (2) 
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Pr = fEr (closure relation). 

For the material quantities, p is density, m is momentum density, p = (7 — l)e is pressure, E 
is energy density, and T is temperature. For the radiation quantities, Er is energy density, 
Fr is flux, Pr is pressure, and f is the variable tensor Eddington factor. In the source terms, 
aa is the absorption cross section, a^ is the scattering cross section, and at = Ca + <Js is the 
total cross section, f is used to close the hierarchy of radiation transport moment equations 
such that: 



Er- 



Idndu, 



nldndu. 



n (X) nldndu. 



The above non-dimensionalization admits the equation of state p = pT [RTao/a^). If one 
assumes physically relevant reference quantities (e.g., gas constant R = 8.31 x 10^ erg K^^, 
Too ~ 10 K, and a^ ^ 3 x 10^ cm s'^), then one finds that C ^ 10^ RT^/al^ ^ 1, and 
p = pT. If one further assumes that poo ~ 10~^ g cm~^, then P ^ lO^^*'. However, the choice 
of these values is arbitrary. These quantities were chosen because they are similar to the 
physical conditions associated with well known test problems (e.g., Su-Olson non-equilibrium 
diffusion as well as subcritical and supercritical radiative shocks). It is important to distin- 
guish that while C and P set the physical scale for the overall system, P is different from the 
ratio Pr/p which can be greater than unity in certain regions of a physical system. 

For the above system of equations, one has assumed that scattering is isotropic and coherent 
in the comoving frame, emission is defined by local thermodynamic equilibrium (LTE), and 
that spectral averages for the cross-sections can be employed (gray approximation). The 



source terms are given by the modified Mihalas-Klein description 31 , 33 which improves 
upon the original Mihalas-Klein description 38 because it maintains important (9(1/C^) 



terms that ensure the correct relaxation rate to thermal equilibrium. See Section 4 of Lowrie 
et al 1999 for a physical description of what these source terms mean. 
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Lastly, because the focus of this paper is exploring the dynamical properties associated with 
resolving a system of stiff hyperbolic balance laws, one has assumed a gray approximation. 
Under this assumption and given the algorithmic machinery built in this paper, a mixed 
frame approach is straightforward to implement. However, this formulation becomes in- 
creasingly complicated for problems defined by multigroup physics where spectral averages 
for the cross-sections cannot be employed. In this context, a comoving frame approach is 
attractive. 



2.1 System of Equations in One Spatial Dimension 

This paper is concerned with formulating the algorithmic ideas for a hybrid Godunov method 
for radiation hydrodynamical problems and develops the method in only one spatial dimen- 
sion. However, the technique should generalize to multidimensions and such work will be 
the subject of a future paper. If on neglects transverse flow, then in one spatial dimension 
the equations can be written as: 
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The quasi-linear form of this system of balance laws is: 
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is specific enthalpy. The Jacobian A has 



eigenvalues: \ = u,u± a, i/^/^C. However, one must account for how the stiff momentum 
and energy source terms affect the hyperbolic structure. 



2.2 Asymptotics 

Before describing the method, it is instructive to consider the properties of Equations [9} 
13 in various limits. Following the asymptotic analysis of Lowrie et al 1999, one considers 
the hyperbolic-parabolic behavior associated with the above system of equations. For non- 
relativistic flows, 1/C = 0{e) where e ^ 1. Assume that there is a moderate amount of 
radiation in the flow and that scattering effects are small where (Js/cTt = C^(e)- 



The optical depth £ is a useful quantity for classifying the limiting behavior of a system that 
is driven by radiation hydrodynamics: 



C 



(JfUX O^il^Cmax -^n 



(14) 



Optically thin regimes are characterized by £ < 0{1), whereas optically thick regimes are 
characterized by £ > 0{1). In optically thin regimes (free streaming limit), radiation and 
hydrodynamics decouple such that the resulting dynamics resemble an advection process. 
In optically thick regimes (weak/strong equilibrium diffusion limit), radiation and hydrody- 
namics are strongly coupled and the resulting dynamics resemble a diffusion process. Given 
these definitions, one makes the following assumption £ = imat/^t = ^mat o"t, where At is the 
total mean free path of the photos and i^^t = C^(l) is the material flow length scale 33 . 



2.2.1 Free Streaming Limit: (JajCrt ~ 0(e) 



In this regime, the right-hand-side of Equation H is negligible such that Equation [9] is 
strictly hyperbolic. Moreover, / — !■ 1 and the Jacobian of the quasi-linear conservation law 
in Equation [12] has eigenvalues X = u,u ± a, ±C. 
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dm d f m^ \ „ ,-,^^ 

f^ + ^ (CF.,) ^ 0. (18) 

^■^' + I- (CBJ = 0. (19) 



9t dx 
2.2.2 Weak Equilibrium Diffusion Limit: (Ta,cTt~C(l) 



One obtains this limit by defining Oa and Ot to be of order unity in Equation 11, matching 



terms of like order, and combining the resulting expressions. From the definition of the 
equilibrium state, Er = T^ + 0{e) and Fr = eu{Ej. + P^) — -^- Therefore, the system 



IS 



parabolic and resembles a diffusion equation, where / — )■ 1/3 

dp d 



dm d f m? \ , , 



dt dx \ p J dx"^ \ 3(Ta 

where: 

p* =p + ¥Pr=p+ -PT^ E* = E + PE^ = E + PT^ e* = e + PT^ (23) 



If one only considers the left-hand-side of Equations [20p2l the Jacobian A^^^ has eigenvalues 
A = M,M ± a*, where a* is the total radiation modified sound speed, a* represents the com- 
bined influence of material and radiation (i.e., effects associated with source terms as well as 
multiscale wave speeds) and is the propagation speed which characterizes the above reduced 
dynamical system. Therefore, this quantity is different than the sound speed associated with 
just the material components of radiation hydrodynamics. 

. *N2 * ^ P*P*e* , _dp* , _ dp* 

(«)=Pp + ^^' Pp = ^' Pe^ = 9^- (24) 

Additionally, the radiation modifies the material such that [33] (z) if the equation of state is 
thermally stable, then a and a* are real-valued quantities, [ii) if the hydrodynamic system 
without radiation is hyperbolic, then the reduced system is also hyperbolic, and (iii) for 
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7-law gases, a* > a. If one assumes a 7-law gas, then one can evaluate a* according to the 



following relationship 40 



(«y^I>^r(p + P/E,)^ ^J^ + 20C + 16C^) ^^^_ p^^ 



^7-1 



120(1 + P 



where T is the Gruneisen coefficient and ^ admits the following limits 33 , 40 60 



^ -)■ =^ r -)■ 7 (no radiation), (26) 

4 
^ — )■ 00 ^ r — ;■ - (pure radiation). (27) 



2.2.3 Strong Equilibrium Diffusion Limit: aa,0't ~ 0{l/e) 



One obtains this limit by defining o"a and at to be much larger than order unity in Equation [TT 
and following the steps outlined for the weak equilibrium diffusion limit. From the definition 
of the equilibrium state, E^. = T^ + 0{e) and F^ = eu{Er + Pj-)- Therefore, the right-hand- 



side of Equations 20 22 is negligible such that the system can be considered hyperbolic. The 



Jacobian of the quasi-linear conservation law has eigenvalues X = u,u±a*, where / — )• 1/3. 

| + |m^0, (28) 

dm d f m? \ , , 

^^•+^f(£-+p-)!^l=0. (30) 



dt dx \ p 

2.2.4 Isothermal Limit 

Lastly, Lowrie et al 1999 evidenced an additional limit which comprises a large portion of 
parameter space: 0{e) < C < 0{l/e). This isothermal limit has some dynamical properties 
in common with the weak equilibrium diffusion limit, but its defining characteristic is that 
the material temperature T{x,t) is constant. Additionally, acoustic waves associated with 
this parameter regime propagate according to the isothermal sound speed (a*)^ = p/p. This 
limiting behavior results from temperature deviations on the slow material fiow scale quickly 



reaching equilibrium on the fast radiation flow scale 33 . This paper will later show how the 



hybrid Godunov method reproduces this isothermal limit. 
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3 Overview of Algorithm and the Nike Code 

In radiation hydrodynamics, there are three important dynamical scales: the speed of sound 
(material flow), speed of light (radiation flow), and speed at which the source terms inter- 
act. When the matter-radiation coupling is weak (free streaming limit), the source terms 
define the slowest scale and the speed of light defines the fastest scale. However, when the 
matter-radiation coupling is strong (equilibrium diffusion limit), the speed of sound defines 
the slowest scale and the source terms as well as the speed of light define the fastest scale. 
Therefore, one must contend with stiffness arising from some waves propagating at the speed 
of light as well as source terms that can have large magnitudes. 

From a mathematical perspective, radiation hydrodynamics can be thought of as a system of 
hyperbolic balance laws with the following dual multiscale behavior: (i) multiscale behavior 
associated with the hyperbolic wave speeds, such that c/ooo ~ 10^, which causes breakdowns 
in monotonicity and (ii) multiscale behavior associated with source term relaxation, such 
that S/ttoo ~ [0, 10^], which causes breakdowns in stability. Despite these dual behaviors 
being different, they both are sources of stiffness that influence the temporal resolution of 
the problem. Rigorous definitions for monotonicity and stability are presented in Matus & 
Lemeshevsky 2009. Given these variations, one desires a numerical technique that treats the 
material flow (slowest hyperbolic scale) explicitly, radiation flow (fastest hyperbolic scale) 
implicitly, and source terms (slow and fast relaxation terms) semi-implicitly. 

The hybrid Godunov method was chosen over two other algorithmic ideas (fully implicit 
methods and PnPm schemes) because of the method's reliance on familiar notions from 
Godunov-type techniques, ability to be easily implemented with minimal computational 
overhead, and accuracy across a wide range of physical behavior. One could have built 
a fully implicit method that advanced time according to the material flow scale, thereby 
circumventing any stiffness related to some hyperbolic waves propagating at the speed of 
light as well as matter-radiation coupling terms. Such an approach was not pursued be- 
cause these methods often have difficulties associated with conditioning, invoke root finding 
techniques when nonlinearities are present, are computationally expensive because of matrix 
inversion, are usually built into central difference schemes rather than higher-order Godunov 
methods, and have trouble sharply resolving discontinuous solutions associated with ad- 
vection and shock phenomena. PnPm schemes are a new family of arbitrary higher-order 
accurate numerical methods for hyperbolic conservation laws. The A^ designates the degree 
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of the polynomials used for test functions in a quasi- WENO (Weighted Essentially Non- 
Oscillatory) technique to spatially reconstruct cell-centered quantities (ll 12 . This WENO- 



like reconstruction is found in ADER (Advection-Diffusion- Reaction) schemes 12 , 55 . The 
M designates the degree of the polynomials used for flux and source term computation in a 
local space-time DG (discontinuous Galerkin) finite element procedure for temporal evolu- 
tion 11,12 . A general feature of the PnPm family of schemes is that they contain classical 
higher-order finite volume methods (A^ = 0) as well as DG methods (M = A^). Due to 
the inherent stiffness of radiation hydrodynamical problems, some authors have suggested 



abandoning classical Godunov methods in favor of DG techniques 31 , 32 . However, such 
methods would be computationally more expensive than the hybrid Godunov method dis- 
cussed here and more difficult to implement in existing codes. Additionally, PnPm schemes 
would not be able to contend with stiffness arising from some waves propagating at the 
speed of light and some modification is required. Nevertheless, PnPm schemes have been 



successfully applied to problems in resistive relativistic MHD 13 and could be applicable 
to other characteristically stiff problems. 

3.1 Effective CFL Condition 

Resolving the evolution of material quantities is the primary interest when solving prob- 
lems in radiation hydrodynamics. The hybrid method computes (p, m, E) to second-order 
accuracy as the entire algorithm is advanced according to an effective CFL condition: 

A< = ^^4^-. (31) 

maxi(|ui| +aeflf,i) 

where At is the time step, u G [0, 1] is the CFL number. Ax = (xmax — a^min)/Accii is the 
grid spacing or spatial resolution for a given number of computational grid cells A'ceii, and 
maxj(|Mj| + aes,i) is the maximum material wave speed over all grid cells. Furthermore, 
Oeff is an estimate of the radiation modified sound speed that is obtained by carrying out an 
effective eigen-analysis of the material Jacobian. This analysis is presented in a later section. 

3.2 Algorithmic Steps 

After defining At, the algorithm loops over the following steps: 

1. Backward Euler Upwinding Scheme - implicitly advances the radiation quantities from 
time in to time tn+i'- (-E", F") — )■ (i?""*"^, -F""*"^) and handles breakdowns in monotonic- 
ity related to the multiscale hyperbolic wave speeds 
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2. Modified Godunov Predictor Scheme - couples stiff source term effects to the hyperbohc 
structure of the balance laws for the material quantities and uses effective piecewise 
linear extrapolation to spatially reconstruct material quantities at the left/right sides 
of cell interfaces i± 1/2: U^mij^i/2 

3. Flux Function - evaluates the passage of material across cell interfaces using left/right 
material states f^^/ij 2+1/2 ^^^1 an approximate Riemann solver 

4. Modified Godunov Corrector Scheme - semi-implicitly advances the material quantities 
from time t„ to time tn+i- (p", m^-, -E") — )■ {p^^^-, m^^^-, E"'^^) and handles breakdowns 
in stability related to multiscale source term relaxation 

5. Apply boundary conditions 

6. Compute next time step At 

In the above expressions, U represents all of the conserved quantities, U^ represents the 
conserved radiation quantities {Er, Fr), and U"^ represents the conserved material quantities 
{pjfn, E). i represents the location of a cell center, i ± 1/2 represents the location of a cell 
interface to the right/left of z, and n represents the time discretization. Details about each 
step are explained in later sections. 

3.3 Numerical Properties 

When solving a system of hyperbolic balance laws with stiff source terms, there are six prop- 
erties that a method must satisfy in order to produce confident numerical solutions. The 
first three properties pertain to solving hyperbolic conservation laws and are: consistency 
(truncation error Ta — )■ as At, Ax — )■ 0), stability (numerical solution remains bounded 
even in the limit of n — > 00 and At — )■ 0, where time t„ = nAt is finite), and accuracy 
(solution to a finite difference equation ip^ converges to the true solution ip at some rate as 
At, Ax — )■ 0). These properties are guaranteed because the hybrid method is based on the 



classical Godunov scheme 12,27,41,57 



If there are source terms, stiffness, and multiple scales associated with the system of equa- 
tions, then a method must also be: coarsely gridded (grid size does not need to resolve 
small scale features attributed to the source terms, therefore At, Ax are based on the asso- 
ciated homogeneous hyperbolic conservation law), well-balanced (method preserves steady 
states), and asymptotic preserving (method gives the correct asymptotic behavior during 
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the relaxation of the system by the stiff source terms) [12]. The first property is satisfied 
because the radiation quantities are imphcitly advanced and the material quantities are semi- 
implicitly advanced according to a time step At which is defined by a conventional CFL con- 
dition that is governed by the material flow scale. The second property pertains to problems 
where there is a balance between flux and source terms and the system exhibits no temporal 
change. This property is satisfied in the modified Godunov corrector scheme when the error 
estimate e(At) = 0, such that | (s"'^, f/'^'«+i) + 5™(f/"^.«, f/'-'"+i)') = V ■ F^^.'^+Vs, Lastly, 
Miniati & Colella 2007 showed their modified Godunov method to be asymptotic preserving 
by looking at the truncation error for a general system of stiff balance laws and examin- 
ing the behavior of the eigenvalues Xcs for the effective Jacobian Aes. The last property is 
satisfied by a^s obtaining the correct asymptotic behavior as one changes the parameters 
{C, P, o"a, at}. This effect is examined in a later section. 

4 Backward Euler Upwinding Scheme 

This section presents the implicit scheme that advances the radiation quantities [Er, Fj.) 
according to the material flow scale. The stability of explicit schemes (e.g., Godunov-type 
methods) is governed by the CFL condition which restricts the allowable time step according 
to the fastest characteristic speed. However, if the overall hyperbolic system consists of mul- 
tiscale wave speeds, then explicit schemes can become inefficient. Radiation hydrodynamics 
(characterized by the speed of light c and material sound speed a^o as well as the speed at 
which source terms interact) is an example of a system that exhibits dual multiscale be- 
havior. This section addresses the multiscale behavior associated with the hyperbolic wave 
speeds. 

4.1 Implicit Schemes and TVD Conditions 

It is well known that numerical difficulties arise when applying implicit schemes to systems 
of hyperbolic conservation laws, particularly when discontinuities are present. Moreover, 
in order to ensure non-oscillatory solutions when using linear implicit higher-order time in- 
tegration methods, one must impose time step restrictions because of the total variation 



diminishing (TVD) condition 14,18,19 . Following the presentation in 20 , the total varia- 



tion of a mesh function v is defined by the following equation: 

oo 

TV{v)= Y, l^m-^.l, (32) 
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where the computational grid is indexed by i. One uses this relation to refer to a numer- 
ical scheme {v"'^^ = S o t>") as being TVD, if TV(f"+^) < TV^(f") where f" denotes an 
approximation to the true solution at time t„. This idea is important to numerical methods 
because Harten 1983 proved that: (z) a monotone scheme is TVD and {ii) a TVD scheme 
is monotonicity preserving. In this context, a numerical method is monotonicity preserving 
if (z) no new local extrema are created within the spatial domain of the numerical solution 
and (ii) the value of a local minimum is non-decreasing which is to say, the value of a local 
maximum is non-increasing. 

From a TVD perspective, Duraisamy & Baeder 2007 presents ratios of the maximum al- 
lowable time step for a few implicit schemes with respect to the maximum allowable time 
step for the explicit forward Euler scheme. From these ratios, one sees that linear im- 
plicit higher-order schemes are impractical when solving flow fields with discontinuities 
such that (At)?^/(At)exp ~ C^(l)- Only the implicit backward Euler method allows for 
(At)^/(At)cxp — !■ cxD [14]. This condition is important when numerically solving radia- 
tion hydrodynamical problems because the backward Euler scheme, which is unconditionally 
TVD stable, can handle the high degree of stiffness associated with the following time scales: 



(At)P^Ep (At) 



c 



,^ , >7T^~ — -10- (33) 

Furthermore, the above statement is consistent with the Godunov order barrier theorem, 
which states: linear numerical schemes for solving partial differential equations that have the 
property of not generating new extrema (i.e., a monotone scheme) can be at most first-order 
accurate ITtI. 



4.2 Linearity and Choice of Algorithm 

If one thinks about these numerical difficulties from a physical perspective, then one can 
intuitively understand why implicit higher-order schemes are not viable approaches for a 
coupled hyperbolic system where some waves move at the speed of sound and other waves 
move at the speed of light. Clearly, there are going to be issues with higher-order schemes 
that do not enforce limiting conditions because the fast waves will have traversed a large 
distance in the computational domain, altered the reconstruction, and caused significant 
oscillations. Therefore, one needs to ensure monotonicity in a higher-order implicit scheme 
that attempts to solve a stiff problem. However, higher-order implicit updates for the ra- 
diation quantities are unlikely because in order to preserve monotonicity one must invoke 
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TVD limiting or artificial viscosity conditions. When applying limiting techniques to the 
implicitly defined linear radiation subsystem , one forces the subsystem to be highly non- 
linear. This approach makes the problem difficult to solve, especially in multiple spatial 
dimensions where convergence is not guaranteed. When using artificial viscosity, one adds 
a competing parabolic/diffusion process. This is a dangerous approach and will likely give 
the incorrect solution in certain parameter regimes because radiation hydrodynamics already 
exhibits parabolic relaxation. Therefore, one is forced to use the first-order backward Euler 
technique. Yet, one way to obtain better resolution while maintaining monotonicity is to 
use adaptive mesh refinement (AMR). Such a technique would be invaluable for accurately 
depicting discontinuous radiative phenomena and is an future area of algorithmic research. 

Since backward Euler-type schemes naturally handle stiff problems and since a primary goal 
of this work is to resolve the material components of radiation hydrodynamics, a backward 
Euler-type scheme provides the suitable framework for defining the implicit update step. 
The backward Euler framework can be cast in the following form: 

jjn+l ^ jjn ^ ^^L (f/"+l) (34) 

where L is a linear operator that defines spatial differences and numerical fluxes. The 
linearity of this operator is important because the radiation subsystem has no nonlinearities 
with respect to Er or Fr. Nonlinearities only arise in the material quantities (p, m, E) which 
are held at time t„. Lastly, Duraisamy & Baeder 2007 describe how to approximate the flux 
integral at cell interfaces which takes the form: 

F^+l/2 = ^ / F (f/(a:,+i/2, t)) dt, (35) 

such that one constructs an implicit flux function. With this idea in mind, an HLLE frame- 
work was chosen for this algorithm because the HLLE flux function is a simple approximate 
Riemann solver with numerical properties that guarantee conservation and good behavior for 
isolated discontinuities as well as smooth solutions. Choosing to implicitly advance the radi- 
ation quantities following the backward Euler HLLE scheme ultimately splits the radiation 
and material components even though the modifled Godunov scheme that semi-implicitly 
advances the material quantities fully maintains the coupling between the material compo- 
nents and source terms. However, this procedural choice does not decrease the accuracy or 
stability of the respective schemes. Moreover, the sum of the total energy (E + Er) across all 
space X is conserved over time. These results are evidenced in the numerical tests presented 
later in Sections 8-10. 
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4.3 HLLE Framework 

The HLLE scheme is based on estimating the minimum and maximum wave speeds (smjn, Smax) 
that arise in the Riemann problem: t/i+1/2 = 'R-{UL,i+i/2,UR,i+i/2)- The numerical flux is 
calculated from: 

i^™(7^([/,, Un)) = Uj'l + :Fn) + \ ( '-■-+_'-^A {T, - T^) , (36) 

Z -^ V ■'max ''inin / 

J^L = F{Ul) - s^JJl, J^r = F{Ur) - s^^^^Ur, (37) 

By combining these relations, one arrives at: 

FH'f/^iniUL, Ur)) = \{{l + C') {F{Ul) - s^JJl) + (1 - C') {F{Ur) - s^^JJr)) . (38) 



where C^ = (s^ax + ■5min)/(smax " Srain)- Inserting C** into Equation 38 results in the more 
familiar relation: 

TTrHLLE/T^/rr TT \\ ^rnax-T [U l) — Srain-T [U RJ + SrninSmax[U R — U LJ , . 

i'i+i/1 [K-Wl, Ur)) = — . (39) 

■'max "Jinin 



However, Equation [38] is the form used by the algorithm's backward Euler update step. 
Defining the left/right states of the Riemann problem according to a first-order accurate 
(piecewise constant) reconstruction, forces the HLLE flux function to become: 

Fl'^^^inUL, Ur)) = Fll^^^inU,, f/,+i)), (40) 

such that the exact integral formulation of the conservative differencing is: 

^n+i = f/f - 1^ (F,+l/2(7^(^7„ f/,+i)) - F,_l/2(7^(f/,_l, U,))) + AtSiU^. (41) 

The choice of using a first-order accurate (piecewise constant) reconstruction was made to 
simplify the calculation and because higher order accuracy is not required since backward 
Euler-type schemes are already first-order accurate. One makes the above explicit HLLE 
scheme implicit by defining the variables in the flux and source terms to be at time tn+i- 



pHLLE 
-^1+1/2 2 



2 ((l + ^r+1/2) {F{Ur') - Sr^JJr') + (1 - Cr+1/2) {F{U?:,') - SmaxC/;VV)) • (43) 
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4.4 Applying the Backward Euler HLLE Scheme 

Consider only the radiation part of the equations for radiation hydrodynamics (i.e., Equations 
^andpl) (49]. This simpler system is termed the radiation subsystem and the variables, fluxes, 
and source terms are: 



u'= : , r\u) = _: , (44) 





s'iu) = I :;. 1 = I L r /. n.n..„;'^^ V .„. 'I y^ h (45) 



c 



-(Tt [Pr 



(l+f)mEr 



CFr \ 




CfEr J ' 


\m f p il+f)mEr\ 

^>pc y^- pc j\ 


)+a^^{T^-Er)_ 





pC 

where the eigenvalues of the radiation subsystem in the free streaming limit (o"a, at ~ 0{e)) 
are A^ = ±/^/^C. Given that the HLLE scheme uses minimum and maximum wave speeds 
to compute fluxes at cell interfaces i + 1/2, one defines the following equations: 

.1/2 _ .1/2 

-Smin = A^ j = — Jj L, Smax = '^R^i+l — Ji+V^^ ^i+l/2 = ~ij2 7172' ^ ' 

Ji+1 + Ji 

where / arises from the closure relation Pr = f'Er and is either a user defined quantity or 
obtained by solving the radiation transport equation. If / varies spatially, then C^_^_^,2 is 
non-zero. Defining or computing f{x,t) precedes the backward Euler update of f/'". For all 
of the numerical tests presented in this paper, / is assumed to be spatially and temporally 
constant, thereby setting C^^^,^ = 0- Future work will update f{x,t) at each iteration by 
solving the radiation transport equation. 

4.5 Matrix Equation for the Radiation Components 



Inputting f/'-.'^+i, F'^([/^''"+i), and ^'■(f/'":"^ [/^•."+i) into Equations |42| and g gives the fol- 
lowing difference equations: 
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e: 



m+l 



-dJi + cu,^ 



1 + dl ( 1 + CiVi/2) /■ /' + rfl (1 - ^7^1/2) /■ /' + d2<Ja 



(47) 



^2 (^g - <ys) (1 + .A) (mf )' 



'M+l 



K 






rfl l + C^l/2 -rfl l-Ctl/2 - 



^2 (o-a - crs)mf 



1/2 






pn+l 
r,i+l 



n+1 



K 



pn+l 



77in+l 



pn+1 



pn+1 

pn+1 
r,i+l 



-d^{l + Cti/2 j /.-I 
'^1(1 + ^^1/2)/.-^? 



dl l + C^l/2 .A-rfl l-Ctl/2 /. 



^20-* (1 + /i) n^? , d2(Tam'2 



Pi 



Pi 



1 + rfi ( 1 + Cf+1/2) /■ /' + rfi (1 - Cti/2) //^' + d20t 



di [^ - C^i+1/2) ft+i 
-di h - Q+1/2J fi+i 



F 



p?<c 



(48) 



where di = AtC/2Aa;, ds = AtC, and Tf = p^/p^ = (7 - 1) 



El 



~'Tp?Y 



If one 



attributes variables 9 and to the bracketed quantities in the above relations, then the 
difference equations can be written as: 



e,E!^;_\ + e,F-t\ + 9,e^;' + 9,f:;' + 9,e-^i, + ^6F;:i\ = e^, + 9,, 



(49) 
(50) 



Since there are no nonlinearities in the radiation quantities for which root finding (e.g., New- 
ton's method) must be implemented, one casts these equations into a sparse matrix format 
that can be solved exactly with basic linear algebra techniques (i.e., Gaussian elimination 
and back substitution). The formulation Ax = b has dimensions dim(A) = 2N x 2N, where 
elements {i?r,i, -^r,i, -E-r.Af, -^r,Af} account for the boundary conditions. Efficient matrix op- 
erations can be performed across the relevant radiation hydrodynamical parameter space 
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because the matrix A is diagonally dominant. Four types of boundary conditions (periodic, 
outflow, reflecting, and inflow) were applied to the numerical tests presented in this paper. 
The ways in which these boundary conditions modify the matrix equations are shown in 
Appendix 1. 

It is important to note that the backward Euler upwinding scheme naturally accounts for 
radiation dominated problems {Pr > p ot E^ > T^). As was discussed earlier and will be 
shown in the following section, the modified Godunov method preserves the isothermal limit. 
When Pr increases with respect to p, material effects as well as updates from the modified 
Godunov scheme become less important. For such physical phenomena, the backward Euler 
upwinding scheme becomes the dominant numerical method in the overall algorithm. If 



one examines Equations |2H6] as well as Equations |44] and |45| then one notices that P is 
only associated with the source terms of the material components. To this end, radiation 
dominated systems are naturally handled by the backward Euler upwinding scheme and is 
evidenced in the numerical tests presented in this paper. 



5 Modified Godunov Predictor Scheme 

Given that the radiation quantities [Ej., Fj.) are at time t„+i, one computes the flux divergence 
(V ■ i^"^)"+i/2 for the material quantities (p, m, E) that are at time t„. Following the analysis 
of Trebatich et al 2005, Miniati & Colella 2007, and Sekora & Stone 2009, one applies 



Duhamel's principle to the quasi-linear system of balance laws in Equations 12 and 13 for 
only the material components. This technique defines the following system that locally 
includes in space-time the effects of the stiff source terms on the hyperbolic structure: 

r)TTrn / f)TJ^ \ 

—^ = l(r]) i -^L-rf- + ^'"(f/'^''^, ^7'^'"+^) J , (51) 

where ^^^ = ^^^ + u^^^ is the total derivative, X is a propagation operator that projects 
the dynamics of the stiff source terms onto the hyperbolic structure, and A^ = A^ — ul 
is the Jacobian for a Lagrangian trajectory with / being the identity matrix. Since the 
predictor scheme is a first-order accurate step in a second-order accurate predictor-corrector 
method, one chooses t] = At/2 and the effective balance law becomes: 

-^ + A'^s-g^ = X(At/2)5'"(f/'"'", f/"'"+^), (52) 

where A"^ = I{At/2)A^ + ul. In order to compute X, one first computes V(7™S'™(f/). 
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5.1 Applying the Modified Godunov Predictor Scheme 



If one only considers the material component in Equations [9 11, then the variables, fluxes, 
and source terms are: 



U'' 




m 



\ 



F"'{U) 



+ p 



(E + p)^ J 



^™(f/) 



Therefore: 



f ° ^ 




/ 




-FS^ 


= 


-P 


\ -FCS^ ) y -PC 


a 



-ot F, 





(l+/)m£;,. 



+ a,^(T^-E, 



o^ij' - E^) + (a, - a.)f F. - (i±{^ 



/ 



Vc/™5'"(f/) 







\ 



V -pc^f -pcs"^ -pc^l / 

where the partial derivatives are: 

-ot (1 + /) mE^ (Tarn (T^ - Er) AaamT^ 






jF 
^m 






Si 






qE 



p2C p2C + pC 



(7-1) 



-E m^ 



p2 p'^ 



pC 



+ 



pC 



+ 



pC 



(7-1) 



— m 



4a,mr3 /l 

(7-1) 



pC 



4cr,T3 (7 - 1 



P. 
E m^\ m{aa- o-J ^r , 2 (a^ - a,) (1 + /) m^E^ 



n2 "*" ^3 



+ 



p- p- J p2C p3C2 

3 ^^, _ , ^ /'^^A , ((Tg - cTs) -^r _ 2 (q-g - q-g) (1 + /) mEr 



p- 



pC 



2(^2 



p2C 



4agT3(7-r 



5.1.1 Simplifying Vc7™5™(f/) 



(53) 



(54) 



(55) 



(56) 
(57) 
(58) 
(59) 
(60) 
(61) 



In its current form, V[/mS'™(f/) in Equation 55 leads to a propagation operator X that is 
difficult to work with algebraically. One should focus on the fact that the material mo- 
mentum source term —FS^ is not the dominant factor defining the stiffness associated with 
the matter-radiation coupling. By inspection, —FS^ < 0{1) even in the strong equilibrium 
diffusion limit, where the contributing terms have the following magnitudes P ~ 0(1/C), 
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(Ta,s,t ~ C(C), Fr < C(l), (1 + f)mEr/pC < C(l), m/pC ~ C(l/C), and (T^ - Er) < C(l). 
Additionally, one finds that the derivative of the material momentum source term with re- 
spect to the conserved variables has the following magnitude — IP5'f „,£} ^ ^(1)- Therefore, 
— P^-^ can be included like a body force (i.e., gravity). 

It is the material energy source term —FCS^ that defines the stiffness associated with the 
problem. By inspection of the contributing terms, — PCS*"^ < (9(C) in the strong equilibrium 
diffusion limit. Additionally, one finds that the derivative of the material energy source term 
with respect to the conserved variables has the following magnitude "PCS'^^^^i < (9(C^). 
Therefore, one only needs to use —FCS^ to define V;7mS'™(f/), such that: 



/ 



Vu^S'^iU) 







:<£' 



\ -FCS^ -PCS, 






-FCS§ 



\ 
J 



(62) 



V(7mS'™'(f/) is further simplified by examining SF ^^^ in the equilibrium diffusion limit and 
neglecting terms that have magnitudes of or less than (9(C). Therefore: 



5f 


—> 


4a,T3 (7 - 


-1) 




-)■ 


4a,T3 (7 - 


-1) 


<.E 
^E 


-)■ 


4a,T-^ (7 - 


-1) 



-E m" 



-m 



(63) 
(64) 
(65) 



It is important to note that these partial derivatives have the same stiff magnitude 4o"aT^ (7 — 1). 
This insight simplifies algebraic manipulation. 

If VumS'^^U) is diagonalizable, then VumS'^{U) = RDR-\ Here, D = diag(0,0, -PCSf) 
and i? is a matrix whose columns are the right eigenvectors. Below, one sees how the stiff 
magnitudes cancel out: 



R 



( 



\ 





5f 








1 






\ 



1/ 



R-^ 



( 



\ 



s§ 



si 



s§ 



s§ 



o\ 



1 



(66) 
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5.1.2 Propagation Operator X 

The propagation operator X is defined in Miniati & Colella 2007. Since one is considering a 
modified Godunov scheme with a predictor step of At/2: 






1 



At/2 



At/2 



^rWumS"^(U)^, 



T 



o\ 











V 



V 



(a-l)|| (a-l)fl a 



J 



1 





1 



o\ 





1 



« 



a] 



a 



J 



(67) 
(68) 

(69) 



where a = (l - exp(-PC^|At/2)) / (PC^|At/2). Since S§ > across all relevant pa- 
rameter space, < a < 1. This property is important when considering stability and the 
subcharacteristic condition which is discussed later in the paper. 



5.2 



Effective Material Jacobian ^4^ 



Before applying X to A™, it is useful to understand that moving-mesh methods can be accom- 
modated in non-relativistic descriptions of radiation hydrodynamics whenever an Eulerian 
frame treatment is employed. These methods do not require transformation to the comoving 
frame 31 . Since the non-dimensionalization is associated with the hydrodynamical scale, 
one can use Wmesh = u from Lagrangian hydrodynamical methods. 



The effects of the stiff source terms on the hyperbolic structure are accounted for by trans- 
forming to a moving-mesh (Lagrangian) frame A™ = A"^ — ul, applying the propagation 
operator X to A^, and transforming back to an Eulerian frame gives the effective material 
Jacobian A^ = lAf + mJ HI]: 



2 ^ 



-(7-3)« 



(7_l„2 _ „H - (1 - a) (^ + lu^)) -(7 - IW + ai/ + (1 - a) (^ + fn^) 





(7-1) 



(70) 
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which has eigenvalues A^|_q_,_| = u — acs,u,u + agg. Here, the effective sound speed a^s 
(i.e., the radiation modified sound speed) is: 

«efr = -^u' + ai^-l)H+il-a)(T + ^u'^ (71) 

= a— + (l-a)T (72) 

P 

= (a(7-l) + l)^, (73) 

P 

where T = p/p because of the non-dimensionahzation in Section 2 admitting the equation of 
state p = pT^RToo/a^). Here, one notices that if, (T + (7 — l)u^ /2) > across all relevant 
parameter space such that the effective sound speed a^s admits the following limits: 

- FCS§ ^O^a^l => alf^^ -^^^^ + (7 - 1)^ = — (adiabatic) (74) 

z p 

-FCS§ ^ -oo^a^O => als^T = - (isothermal). (75) 



When examining Equations 74 and 75, one sees that the subcharacteristic condition for 



material wave speeds is satisfied, such that 41 



A!? = tx - aad < AS,- =u-a,s< X^" = XTs,o = u< A^ff^^ = u + a,s < X"-' = u + a^d- (76) 

This condition is necessary for the stability of the system and guarantees that the numerical 
solution tends to the solution of the equilibrium equation as the relaxation time tends to 
zero. Additionally, the structure of the equations remains consistent with respect to clas- 
sical Godunov methods. Therefore, the CFL condition max(|A^^|)^ < 1 applies, where 

* = {-,o,+}. 

It is important to note that a^s is different from a* which was defined earlier as the propa- 
gation speed associated with a reduced dynamical system that captures the combined influ- 
ence of material and radiation (i.e., effects associated with source terms as well as multiscale 
wave speeds), aes defines the effective material sound speed that estimates the influence 
of the source terms on the hyperbolic structure of the split material subsystem. Although 
the modified Godunov method is an unsplit finite volume technique that directly couples 
advection-reaction processes, there is an obvious splitting between the material (modified 
Godunov scheme) and radiation (backward Euler upwinding scheme) because of the mul- 
tiscale nature of the hyperbolic wave speeds, agg is designed to handle the former piece 
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of this splitting. In the dispersion analysis of Lowrie et al 1999, a* is plotted for various 
parameters. As the optical depth is increased, a* takes values equal to the adiabatic sound 
speed and then values equal to the isothermal sound speed. At very large optical depths, a* 
takes values exceeding the isothermal sound speed. Because aes is only defined for the split 
material subsystem, not the global radiation hydrodynamical system, aes should only take 
values ranging between the adiabatic and isothermal sound speeds which is what is seen in 



Equations 74 and 75 In situations where large optical depths are involved, radiation is the 



dominant species governing the dynamical behavior and is evolved using the backward Euler 
upwinding scheme for the split radiation subsystem. 

Lastly, the right material eigenvectors R^ (stored as columns) and left material eigenvectors 
L^ (stored as rows) are shown below. Clearly, there are many structural similarities between 
the effective eigen-quantities A^, A^, K^, and L^ and those from adiabatic/isothermal 
hydrodynamics. 
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5.3 Computing Left /Right States 

In the modified Godunov predictor scheme one uses effective piecewise linear extrapolation 
to spatially reconstruct material quantities at the left/right sides of cell interfaces. This 
technique is given by the following relations: 



U 



m,n+l/2 
L,i+l/2 
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u 



m,n+l/2 
R,i+l/2 



a 



i+l 



+ 
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^( Y'^"^^' 
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(79) 



(80) 
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where / is the identity matrix, A^ is the effective material Jacobian, Pa is a slope limiting 
function used to eliminate spurious oscillations, X is the propagation operator, and S"^ is 
the vector of source terms influencing the material equations. Slope limiting is performed 
for each of the material quantities. Although any traditional slope limiter can be used, the 
Nike code employs the extremum-preserving (8,48 as well as the traditional van Leer lim- 



iter (also referred to as the MUSCL limiter). These techniques can be implemented either 
componentwise or across characteristic fields and are illustrated in Appendix 2. 



After reconstructing the material quantities, an approximate Riemann solver evaluates the 
passage of material across each cell interface using left /right states U^j'^^j^^,^ 
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The 



Nike code employs three such flux functions (Lax-Friedrichs, Lax-Wendroff, and HLLE), 
which are presented in Appendix 3. It is important to emphasize that these flux functions 
do not directly account for the influence of radiation on the material quantities. Rather, 
the radiation effects are included via the source terms, propagation operator, and effective 
material Jacobian. 

6 Modified Godunov Corrector Scheme 

The time discretization for the source term is a single-step, second-order accurate scheme 



based on the ideas from Dutt et al 2000, Minion 2003, and Miniati & Colella 2007 [15||41||42 
Given the material system of balance laws, one aims for a scheme that has an explicit ap- 
proach for the flux divergence term V ■ F"^ and an implicit approach for the stiff source term 
5™(f/). 

At each grid point, one solves the following collection of ordinary differential equations: 

S'^iU) - (V ■ F'")"+i/^ (81) 



dt 



where the time-centered flux divergence term is inputted from the predictor step and taken 
to be constant valued. Using Picard iteration and the method of deferred corrections, an 
initial guess for the solution to the collection of ordinary differential equations is: 

f/ = f/™." + At(/ - AtVc7™5'"(f/)|c7™.n,t/'-."+i)"^('5''"(^'"''', f^''""^^) - (V ■ F")"+i/2)^ (g2) 
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where Vu"^S"^{U) has the same functional form as that which was used to define the prop- 
agation operator X in a previous section. Therefore: 

1 

{I-AtVu^S"'{U))=\ 1 |. (83) 

AtPCSf AtFCS^ l + At¥CS§ 



By inverting the above matrix, one finds: 
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1 








{I - AtVu^s{u)y^ = 






-At¥CSf 


1 

-AtPC5^ 




1 



(84) 

y l+AtPCSf l+AiPC5| l+AtPC5| / 

The error estimate e is the difference between the solution obtained from one iteration of the 
Picard technique where U is used as the starting value and the initial guess U: 

e{At) = U""'"" + — (s'^iU, f/"'"+i) + 5"(f/™'", [7"'"+^)) - At(V ■ F™)"+i/2 _ fj_ (§5) 

Following Miniati & Colella 2007, the correction to the initial guess is given by: 

5(At) = (^/- AtVt/-5"(f/)lc>,t/'-."+i) 'e(At). (86) 

Therefore, the material quantities at time t„+i are: 

f/™'"+i = f/ + 5(At). (87) 

Appendix 4 explains how to apply various types of boundary conditions in the modified 
Godunov corrector scheme. 

7 Numerical Tests of Nike 

The hybrid Godunov method was implemented in a new radiation hydrodynamical code 
called Nike. A suite of ten numerical tests were conducted to gauge the temporal and spa- 
tial accuracy of the hybrid Godunov method. Tests were carried out using the same code 
with no special adjustments made for any test and span a wide range of mathematical and 
physical behavior. These tests are grouped into three categories (Hydrodynamics, Radiation 
Subsystem, and Full Radiation Hydrodynamics) which define the dominant physical system. 
The numerical solution is compared with the analytic solution where possible. Otherwise, a 
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self-similar comparison is made. 

The following definitions for the n-norms and convergence rates are used throughout this 
paper. Given the numerical solution q^' at resolution r and the analytic solution u, the error 
at a given point i is e[ = g[ — u. Likewise, given the numerical solution q^ at resolution r 
and the numerical solution g^'"*"^ at the next finer resolution r + 1 (spatially averaged onto 
the coarser grid), the self-similar error at a given point z is e[ = g[ — q^'^'^- The 1-norm and 
max- norm are: 

Li = V'|e[|Ax^ Lmax = max|e[|. (88) 

i 

The convergence rate is measured using Richardson extrapolation: 

_ ln(L.(e-)/L46-+^)) 

ln(Ax7Ax-+i) ■ ^ ' 

For all of the tests presented in this paper, the material fluxes were computed with a compo- 
nentwise van Leer limiter and an HLLE flux function using second-order (piecewise linear) 
spatial reconstruction. Moreover, 7 = 5/3 and the CFL number is i/ = 0.5. This value for 
the CFL number was chosen for convenience but any z/ < 1 is valid. 

8 Hydrodynamics 

In order to determine how Nike's modifled Godunov scheme resolves purely hydrodynamical 
phenomena, the following set of parameters were chosen to force the radiation hydrody- 
namical system towards the free streaming limit. Therefore, Equations [9 11 transform into 



Equations 15 17, The following tests examine processes associated with the Euler Equations. 



Parameters: 



C = 10^ P = 10-20, aa,at = lO-^o, / = 1, eO,F° = lO'^o, 



///\t 

^^ = TT^, V' ^-- = 0' ^-ax = 1, N,,n = [32, 64, 128, 256]. 

When conducting purely hydrodynamical tests on the modified Godunov scheme, one expects 
R = 2.0 for smooth initial data (e.g., Gaussian pulse) and R ~ 0.67 for discontinuous initial 
data (e.g., square pulse). This claim is true for all second-order spatially accurate numerical 
methods when applied to an advection-type problem {ut + aux = 0) [27]. It is important to 
present these hydrodynamical tests because some radiation hydrodynamical codes perform 
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Figure 1: Gaussian pulse in the free 
streaming limit (hydro variables), t = 1. 
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Figure 2: Square pulse in the free stream- 
ing limit (hydro variables), t = 1. 



poorly in this free streaming limit. The results herein demonstrate the hybrid Godunov 
method to be robust and applicable across different parameter regimes. 

8.1 Linear Advection 

For the linear advection tests, the initial values of the conserved quantities are given by: 



/ P' \ 



U{x,0) = 

Parameters: 

IC for Gaussian Pulse: 







m 



( p' \ 



\ 



P 

E^ 

r 



(90) 



/ 



p° = 1 + exp {-{v{x - /i))2) , t; = 20, /i = 0.5, 
M° = 1, p° = 1, tend = 1 (one crossing time) 



IC for Square Pulse: { "}'' < ^ < »■«' = '• ""/^ < ''•'■" < ^' = '■'• 

u" = 1, p^ = 1, tend = 1 (one crossiug time) 

periodic BC. 

As the spatial resolution increases, the numerical solutions in Figures [l] and [2] converge to 
Gaussian and square pulses. Furthermore, Tables 1 and 2 show that the convergence rates 
approach 2.0 and 0.67 for smooth and discontinuous initial data. 
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Table 1: Convergence analysis for a Gaussian pulse in the free streaming limit (hydro 
variables). Errors were determined by comparing the numerical/analytic solutions, t = 1. 
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Table 2: Convergence analysis for a square pulse in the free streaming limit (hydro 
variables). Errors were determined by comparing the numerical/analytic solutions, t = 1. 
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Table 3: Convergence analysis for the hydrodynamical linear mode {u — a). Errors were 
determined by comparing the numerical/analytic solutions, t = 1. 



8.2 Propagation of Linear Modes 

One of the most discriminating hydrodynamical tests is the propagation of linear modes on 
a periodic domain. One initializes eigenfunctions of sound and contact waves in a uniform 
medium by using the right material eigenvectors R^ and setting the wavelength of each 
small amplitude perturbation A equal to the size of the domain (L = Xr, 



Xr, 



U{x,0) 



f p^ + A sm{kx) R';;/'''\p^,u^,p^) \ 



E^ 

r 



(91) 



Parameters: 



P 



1, u"" 



(sound waves: m ± a), vP = 1 (contact waves: m), p° = I/7, 



^ = 10 , k = 27r/L, j = {1,2,3}, tend = 1 (one crossing time), periodic BC. 

This choice of parameters causes the material sound speed to take the following value Ogg = 1. 
As seen in Tables 3 and 5, the errors associated with the sound waves {u ± a) are identical 
because the initialization is the same except that the waves propagate in opposite directions. 
As evidenced in Tables 3-5, all of the convergence rates approach 2.0. 



8.3 Sod Shock Tube 

The Sod shock tube is an excellent test of an algorithm's ability to resolve nonlinear behav- 
ior. This well known test consists of two constant states separated by a discontinuity (i.e.. 
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Table 4: Convergence analysis for the hydrodynamical linear mode (u). Errors were 
determined by comparing the numerical/analytic solutions, t = 1. 
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Table 5: Convergence analysis for the hydrodynamical linear mode {u + a). Errors were 
determined by comparing the numerical/analytic solutions, t = 1. 
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Figure 3: Sod shock tube, t = 0.2, iVceii = 256. 



a Riemann problem). 



Parameters: 



U{X < Xint,0) 



/ 1 


\ 







1 

(7-1) 




\ F^ 


/ 



U{xint < a;,0) 



/ 0.125 \ 


0.1 

(7-1) 

T?0 



•^int 



0.5, tend = 0.2, N, 



cell 



256, outflow BC. 



As the Euler equations evolve in time, a rarefaction fan, contact discontinuity, and shock 



wave form 51 , 56 . Figure 3 shows a higher-order accurate solution to the Sod shock tube 



problem, where discontinuities are separated by only one grid cell. 



9 Radiation Subsystem 

In order to develop a reliable and robust numerical method for full radiation hydrodynamics, 
Sekora & Stone 2009 examined the radiation subsystem, which is a simpler set of equations 
that minimizes complexity while maintaining the rich hyperbolic-parabolic behavior associ- 
ated with balance laws that have stiff source terms. If one initially sets the material flow to 



be stationary such that m — )■ in Equations 44 and 45, then the variables, fluxes, and source 
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terms for the radiation subsystem are given by: 

^^Cf ^C.„r-K.). (92) 

^^C/f^-C.F. (93) 

In order to determine how Nike's backward Euler upwinding scheme resolves the radiation 
subsystem, the following set of parameters were chosen. 

Parameters: 

C = 10^ P = 10-20, pO, m°, E^ = 10-20, 

At = ^^, iV,eii = [32, 64, 128, 256]. 

When conducting tests on the backward Euler upwinding scheme, one expects i? = 1.0 for 
smooth initial data (e.g., Gaussian pulse) and R ~ 0.5 for discontinuous initial data (e.g., 
square pulse). This claim is true for all first-order spatially accurate numerical methods 
when applied to an advection-type problem {ut + au^ = 0) l27|. By no means was the 
material flow held stationary {u = 0) throughout the evolution of the following numerical 
tests. All problems were solved with the full dynamical Nike code. The fact that the material 
quantities remain zero in these tests demonstrates the robustness of the overall method. The 
following problems were initialized according to: 

t/(x,0) = (pO,mO,EO,EO,FO). (94) 

9.1 Linear Advection - Free Steaming Limit 



In the free streaming limit, the radiation subsystem reduces to Equations 18 and 19 If one 
takes an additional temporal and spatial partial derivative of the radiation subsystem in 
the free streaming limit and subtracts the resulting equations, then one finds two decoupled 
wave equations that have the following analytic solutions: 

Er{x,t) = Eo{x-f/'^Ct), (95) 

F,(x,t) = Fo(x-f/2o). (96) 

Parameters: 

aa, (Jt = 10-20, / = 1, a;min = 0, Xmax = 1, periodic BC, 
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Table 6: Convergence analysis for a Gaussian pulse in the free streaming limit (radiation 
subsystem). Errors were determined by comparing the numerical/analytic solutions, t = 1. 
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Table 7: Convergence analysis for a square pulse in the free streaming limit (radiation 
subsystem). Errors were determined by comparing the numerical/analytic solutions, t = 1. 



IC for Gaussian Pulse: E°, F° 



exp (— (t>(a; — yu))^) , v = 20, fi 
IC for Square Pulse: E°, F^ = 



0.5, 



1 0.4 < x < 0.6 
otherwise 



Graphically, these initial conditions give rise to numerical solutions with profiles that are 
identical to those shown in Figures [T] and |2} In Tables 6 and 7, one observes first-order 
accurate convergence rates. What is most interesting is how long it takes (i.e., how high of 
a spatial resolution is necessary) to see R —> 1 for a Gaussian pulse in the free streaming 
limit. Yet, this result is not wholly unexpected. In this test, one is examining advective 
behavior. It is widely known that implicit algorithms perform poorly when trying to capture 
the propagation of individual wave modes. The fact that the backward Euler upwinding 
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scheme obtains the correct solution in a stable and accurate manner is a testament to the 
robustness of the overall method. Moreover, resolving individual radiative waves associated 
with the speed of light is not the primary concern of the hybrid Godunov method. Therefore, 
it is of little consequence that such high spatial resolution is necessary to see a first-order 
accurate convergence rate for radiation quantities in the free streaming limit. 

9.2 Exponential Growth/Decay to Thermal Equilibrium 

This test examines the temporal accuracy of how radiation variables are updated in the 
backward Euler upwinding scheme. Given the radiation subsystem and the following initial 
conditions: 

E^. = constant across space, F^ = 0, T = constant across space and time, 

-Fr — )■ for all time. Therefore, the radiation subsystem reduces to the following ordinary 
differential equation: 

^ = Caa{T'-E,), (97) 

which has the following analytic solution: 

E, = T^ + (E°-T^)exp(-C(Tat). (98) 

For E^ < T^ and F^ = 0, one expects exponential growth in Er until thermal equilibrium 
[Er = T^) is reached. For E^ > T^ and F^ = 0, one expects exponential decay in Er until 
thermal equilibrium is reached. Since E^. is spatially uniform for this system, this numerical 
test allows one to examine the temporal order of accuracy of the backward Euler upwinding 
scheme as a stiff ODE integrator. 

Parameters: 

IC for Growth: E° = 1, F° = 0, T = 10, 

IC for Decay: E° = 10^ F° = 0, T = 1. 

From Figure |4| one sees that the numerical solution corresponds with the analytic solution. 
In Table 8, ones sees that at lower resolution the errors and convergence rates are identical 
for growth and decay. However, as the resolution increases one sees disparity between these 
values, such that the growth tests exhibit convergence rates that are less than first-order while 
the decay tests exhibit convergence rates that are greater than first-order. This asymmetry 
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Figure 4: Exponential growth/decay of Er to thermal equilibrium. A'^ceii = 256. 
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Table 8: Convergence analysis for exponential growth/decay in Ej. to thermal equilibrium. 
Errors were determined by comparing the numerical/analytic solutions, t = 10^^ = l/o"aC. 
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is expected for this situation because it is associated with the unconditional stabihty of a 
backward Euler-type integrator. Whenever a problem is defined by dynamical exponential 
growth, an algorithm works hard to maintain numerical stability. However, when a problem 
is defined by dynamical exponential decay, the nature of the problem already ensures stability 
and the algorithm doesn't have to work as hard. Nonetheless, one finds that the method is 
well behaved and obtains the correct solution with at least first-order accuracy for even stiff 
values of the e folding time (y^c — ^)' This result credits the flexibility of the backward 
Euler upwinding scheme as being a robust temporal integrator. 

9.3 Weak Equilibrium Diffusion Limit 

In the weak equilibrium diffusion limit, the radiation subsystem reduces to the following set 
of equations: 

^ = ^^ (99) 

K - -f f^. (100) 

Sat ox 
The optical depth suggests the range of opacities for which diffusion is observed for this 
subsystem. If £ = cij ^diff > 1, then one expects diffusive behavior for at > l/^difr- Addition- 



ally, Equations |99| and |100| set the time scale tdifr and length scale ^difr for diffusion, where 
tdiff ~ ^diff/-^ ^^"^ -^ ~ f'^/'^t- Given a diffusion problem with a Gaussian pulse defined over 
the entire real line {ut — Du^x = 0), the analytic solution is given by the method of Green's 
functions: 

u(., *) = £ /(.)G(., , X, 0)£ = p^J^^exp ( "i^tTD ■ (1°1) 

Parameters: 

<7a) CT* = 40, / = 1/3, Xmin = —5, Xmax = 5, 

A^ceii = [320, 640, 1280, 2560], Outflow BC, 

E° = exp {-{v{x - ^i)f) , t; = 20, /i = 0, 
IC for Gaussian Pulse: I F" = -^^ = 2/^'("-^) ^o 

^ at ax at ^^ 

T^ = Er 

tend = [0, 1, 4, 16] X 10-^ 

One's intuition about diffusive processes is based on considering an infinite domain. So to 
minimize boundary effects in the numerical calculation, the computational domain and num- 
ber of grid cells were expanded by a factor of 10. In Figures |5] and [6| one observes the diffusive 
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Figure 5: Erin weak diflF limit (rad subsys- 
tem), t = [0, 1,4, 16] X 10-^ iVceii = 2560. 
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Figure 6: Fr in weak diff limit (rad subsys- 
tem), t = [0, 1,4, 16] X 10-^ iVceii = 2560. 
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Table 9: Convergence analysis for Er, Fr in the weak equilibrium diffusion limit. Errors 
were determined by a self-similar comparison of the numerical solutions, t = 4 x 10~®. 
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behavior that is expected for this parameter regime. Additionally, the numerical solution 
compares well with the analytic solution for a diffusion process defined over the entire real 



line (Equation 101). However, diffusive behavior is only a first-order approximation to more 
complicated hyperbolic-parabolic dynamics taking place in radiation hydrodynamics as well 
as the radiation subsystem. Therefore, one compares the numerical solution self-similarly 
and sees the expected first-order convergence in Table 9. 

In order to examine self-similar convergence in the weak equilibrium diffusion limit for the 
above tests, the backward Euler upwinding scheme uses parameter values and a fine spatial 
resolution Ax to resolve the photon mean free path At. Therefore, the cell-optical depth £a 
is always fairly small, such that £a = Ax/Aj = Ax at = [1.25,0.625,0.3125,0.15625] for 
(Tt = 40 and A^ccU = [320,640, 1280,2560], respectively. Consequently, the photon mean free 
path is moderately resolved which suggests that the backward Euler upwinding scheme is nu- 
merically consistent but not necessarily asymptotic preserving. This idea was first developed 
in the context of radiation hydrodynamics by Larsen et al 1987 although these authors did 
not use the term asymptotic preserving. In order to explore this idea, consider the following 
hypothetical question - can one say that a conventional Godunov method is asymptotic pre- 
serving for the the system defined by Equations [9 11? If one evolves the system according 



to At ~ 1/C, then one expects diffusion behavior for certain parameter regimes because for 
high spatial and temporal resolution consistency implies asymptotic preservation. However, 
because one has to use such fine spatial and temporal resolution, a conventional Godunov 
method is not asymptotic preserving. 

To demonstrate that the backward Euler upwinding scheme has some asymptotic preserving 
properties, one can use the problem settings that were used earlier but rescale the spatial 
domain, initial conditions, and time for how long the computation runs so that £a increases 
by some significant amount, say 1000. According to this scaling, which is defined by the 
relation tdifr ~ ^dml-^i ^^^ expects the profiles of E^ and Fj, to match the profiles shown 
in Figures [5] and [6] except that the spatial scale will have increased by a factor of 1000. 
Additionally, one expects the magnitude of F^ to decrease by a factor of 1000 because in 



the weak equilibrium diffusion limit, E^ behaves according to Equation |101| and Fr behaves 
according to the following relation: 

FJx,t) = - — —^ = ^——^Er. 102 

^ ' ^ at dx at ADtv^ + 1 ^ ' 
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Figure 7: E^ in weak diff limit (rad subsys- 
tem) defined by large cell-optical depths. 
Red circles designate A^ceii = 320 (£a = 
1250) and black lines designate A^ceii = 
2560 (£a = 156.25). t = [0, 1, 4, 16]. 



o 

X 



0.15 

0.1 

0.05 



-0.05 

-0.1 

-0.15 




-500 -250 



250 



500 



Figure 8: F^ in weak diff limit (rad subsys- 
tem) defined by large cell-optical depths. 
Red circles designate A^ceii = 320 (£a = 
1250) and black lines designate A'^ceii = 
2560 (£a = 156.25). t = [0, 1, 4, 16]. 



Parameters: 



aa, at = 40, / = 1/3, x^i^ = -5000, x^ax = 5000, 
N^eii = [320, 640, 1280, 2560], Outflow BC, 







exp (— (f (x — n)y 



IC for Gaussian Pulse: 



e: 



r at dx 



0.02, /i = 0. 



2 f v' (x- fj.) jj^o 

at f"' 



T-l4 

^end 



Er 

[0, 1, 4, 16]. 

Results for these parameters are shown in Figures [7] and |8| Again, one compares the nu- 
merical solutions self-similarly across different spatial resolutions at t = 4 in Table 10 to 
examine the convergence properties of the backward Euler upwinding scheme. From the 
values seen in Table 10, one initially sees first-order convergence. However, as the spatial 
resolution is refined the error norms stop decreasing. This behavior suggests that for large 
cell-optical depths, the numerical solution has converged on some asymptotic profile where 
higher spatial resolution does not enable further convergence because of a competing process 
associated with the implicit discretization of the radiation subsystem in the backward Euler 
upwinding scheme. Furthermore, if one examines the discretization of the backward Euler 



A Hybrid Godunov Method for Radiation Hydrodynamics 



41 



A^cell 


Li{Er) 


Lo.{Er) 


Li{Fr) 


Lo.{Fr) 


320 


3.9E-2 


l.OE-1 


3.4E-6 


1.3E-5 


640 


2.0E-2 


5.2E-2 


1.8E-6 


6.5E-6 


1280 


8.8E-2 


6.6E-1 


1.7E-5 


1.3E-4 


2560 


8.8E-2 


6.6E-1 


1.7E-5 


1.3E-4 



Table 10: Convergence analysis for Er, F^ in the weak equilibrium diffusion limit defined 
by large cell-optical depths. Errors were determined by a self-similar comparison of the 

numerical solutions, t = 4. 



HLLE scheme, one finds the following difference equation for the radiation energy: 

- (rfi/V2) E^^\ + (1 + 2d,f'' + d^aa) E^^' - {dj'/') El% + ...=El^ + ... (103) 

This difference equation is important because it has the same functional discretization as a 
backward-time centered-space scheme for a parabolic/diffusion equation. Therefore, one sees 
consistency associated with the truncation error and discretization operator. All of the above 
insights suggest that the backward Euler upwinding scheme can preserve certain asymptotic 
limits. However, this does not suggest that the overall algorithm (i.e., hybrid Godunov 
method) is globally asymptotic preserving. Further discussion of this topic is reserved for 
the conclusions at the end of the paper. 



9.4 Strong Equilibrium Diffusion Limit 



In the strong equilibrium diffusion limit, the radiation subsystem reduces to Equations 104 



and 105, which suggest that F^. — )■ for all time and space while E,. = E^. 

dEr 



dt 

Fr = 0. 



0, 



Parameters: 



O'ay C"t — J-U , J — -L/j) Xjjiiji — O, Xmax — t), 

iVceii=[320, 640, 1280, 2560], Outflow BC 



IC for Gaussian Pulse: 



E^ = exp {-{v{x - /i))2) , t; = 20, /x = 0, 

r ( 



at dx 



E'! 



(104) 
(105) 
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Ncoll 


L,{Er) 


R 


LooiEr) 


R 


Ll{Fr) 


R 


Lo.{Fr) 


R 


320 


l.lE-1 


- 


8.1E-1 


- 


9.4E-7 


- 


7.5E-6 


- 


640 


6.1E-2 


0.9 


5.0E-1 


0.7 


6.3E-7 


0.6 


7.0E-6 


0.1 


1280 


3.1E-2 


1.0 


2.6E-1 


0.9 


3.5E-7 


0.8 


4.5E-6 


0.6 


2560 


1.6E-2 


1.0 


1.3E-1 


1.0 


1.9E-7 


0.9 


2.4E-6 


0.9 



Table 11: Convergence analysis for E^, F^ in the strong equilibrium diffusion limit. Errors 
were determined by a self-similar comparison of the numerical solutions, t = 4 x 10~®. 



In this test, the numerical solution remains fixed at the initial distribution because o"a and 
at are so large. However, if one fixed i^is and scaled time such that tdifr ~ ^diff/-^ — 
i^^f^at/fC, then one would observe behavior similar to Figures [s] and [6] This test illustrates 
the robustness of the backward Euler upwinding scheme to handle very stiff source terms. 

9.5 Su-Olson Non-Equilibrium Diffusion 

Due to the complexity of the equations for radiation hydrodynamics, reference problems 
with analytic solutions are rare outside of the free streaming limit. However, Marshak 1958 
presented a time-dependent radiative transfer problem. The test which is described below is 
a variation of the problem presented by Marshak 1958 and was analyzed by Pomraning 1979 
but more extensively by Su & Olson 1996. Therefore, this test is referred to as the Su-Olson 
problem. 

The following test investigates non-equilibrium radiation diffusion as well as the physics 
of matter-radiation coupling and consists of a purely absorbing, homogeneous medium that 
occupies the one- dimensional, semi-infinite space < x < oo. Initially, the medium is defined 
by zero radiation energy Er and material temperature T. Upon beginning the calculation, a 
time independent radiative flux F^"^^ impinges upon the left boundary at x = p4|. If one 
ignores hydrodynamic motion as well as heat conduction and further assumes a diffusion-type 



approximation, then one can write the following simplifled expressions 46 , 54 
dEr{x,t) d fCfdEr{x,t) 



dt dx 



(j;^^^) = ^^a {T\x,t) - Erix,t)) , (106) 



c,{T)^^^^ = -Ca, {T\x,t)-E,{x,t)) , (107) 
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where Cy is the specific heat capacity of the material medium and T is related to the material 
internal energy e = pc^T. If one chooses the following functional form for the specific heat 
Cy = (T^, then Equations 106 and 107 can be written as linear ODEs in Ej. and T^. Here, 



( is an arbitrary constant that if chosen appropriately simplifies how material and radiation 



quantities evolve 46| 54 . Lastly, the boundary and initial conditions are defined by: 

Er{00,t) =0, 



BC: Er{0,t) - -^^^ = AF-^ 



at dx 
IC: Er.{x,0) =T{x,0) = 0. 



(108) 
(109) 



If one sets C = 4 and defines the radiation flux as a dependent variable, then the above 
system can be cast into a form that resembles the radiation subsystem where the material 
temperature temporally evolves with the radiation quantities: 



+ ^ = Caa (T^ - Er) 



dt 



dx 






CfdEr 
at dx 

eCa, [T^ - Er) 



IC: EriO,t)+2FriO,t) =4:F, 



mc 
r 



(110) 

(111) 

(112) 
(113) 



This problem setup along with specific initial and boundary conditions as well as a unique 
choice of parameters causes the radiation subsystem to exhibit well defined dynamics char- 
acterized by the diffusion limit and non-equilibrium behavior. Su & Olson 1996 arrived at a 
semi-analytic solution to the above system by defining dimensionless independent variables 
{X, T) and dependent variables (W, V) such that: 



X 



cr, 



av^; 






T 



eCaJ, 



U{X,r) = Er{x,t)/4Fr, V{X,r) = T\x,t)/4F^ 



mc 
r 



Using these dimensionless variables. Equations 106 109 become |46y54| : 



dT dX^ 
dV 



dT 



W-V, 



BC: W(0,r)--|^^^^^ = l, W(oo,r) = 0, 



(114) 
(115) 

(116) 

(117) 
(118) 
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IC: ^(A'.O) = V{X,0) = 0. (119) 



The first boundary condition in Equation |118 enforces tlie constraint of constant flux on 



tlie left side of the computational domain. This special boundary condition can be imple- 
mented by imposing the time- varying Dirichlet condition U{0, T) pGl. This quantity, which 
is evaluated at A" = 0, is given by l54j: 



uix,T) - 1 - ^ / di.M-re) I ^^^^^^P^=|P) (120) 

sin(A'T2(0 + e2(e)) 




e(l + ee)v/3T4Tl(0 




^i(e) = ^\h+Y^2^ ^2(0 = W(l -0 ( e+ - ), (121) 



e„(0-os-^^^^^^, n^l,2. (122) 

Since the backward Euler upwinding scheme is a first-order accurate method, it suffices 
to use the Trapezoidal Rule Tint(/) = (/(&) + f{ci)){h — a)/2 to numerically evaluate the 



integrals in Equation 120 Higher-order integration techniques were also attempted, but 



these methods had little effect on the accuracy of the solution computed by the backward 



Euler upwinding scheme. One additional difficulty associated with evaluating 120 is that the 
integrands are highly oscillatory for T ^ 1. Therefore, one uses the following asymptotic 
expression whenever T < 10^'^ 46,54 : 



The last element in setting up this problem is updating the material temperature T. In 
previous tests, T was computed by {i) using other material quantities: T = p/p, {ii) defining 
a fixed background temperature: T{x,t) = T°(x,0), or {in) assuming the system to be 
in thermal equilibrium: T^ = Er- Due to the non-equilibrium nature of this problem, 
one implicitly advances T according to the following relation prior to solving the radiation 



variables in the backward Euler upwinding scheme 54 



The following parameters were chosen because of earlier work by by Sekora & Stone 2009, Su 
& Olson 1996, and Reynolds et al 2009. By setting F^'^'^ = 1/4, one normalizes the problem 
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Figure 9: Su-Olson problem with e = 0.1, 
where the sets of curves from bottom to top 
correspond to t = [2.5,7.5,25,75,250] x 
10"^, respectively. A^ceii = 512. 



0.015 



■■■ 32 cells 
— 64 cells 
cells 
cells 
cells 




0.035 



0.055 



0.075 



Figure 10: Convergence results for the Su- 
Olson problem with e = 0.1 and t = 2.5 x 
10^^. The bottom and top sets of curves 
correspond to T"^ and Er, respectively. 



such that Er{x, t) and T^(x, t) directly correspond to U{X , T) and V{X , T) under the appro- 
priate change of spatial and temporal coordinates [54] . Since the Su-Olson problem assumes 
a diffusion approximation, one chooses (JaiCrt = 40 like in the test problem related to the 
weak equilibrium diffusion limit. Reynolds et al 2009 examine the spatial range < ^Y < 35 
with a resolution of 2048 grid cells. These values and the relation X = aavSx define the 
spatial aspect of the computational domain. Furthermore, Su & Olson 1996 tabulate semi- 
analytic solutions for the dimensionless time T = [1,3,10,30, 100]. Therefore, these values 
and the relation T = eCo"at define the temporal aspect of the computational domain. Again, 
this problem was solved using the full dynamical Nike node as calculations were performed 
using e = 0.1, 1.0. 



Parameters: 

a„,a, = 40, / = l/3, F;- = 1/4, 

Xmin = 0, x^ax = 0.5, A^^eu = [32, 64, 128, 256, 512], 

^ , , f e = 0.1 ^t= [2.5,7.5,25,75,250] X 10-6 

r= 1,3,10,30, 100 : <^ [ , , , , J 

[ e = 1.0 ^ t = [2.5, 7.5, 25, 75, 250] X 10-^ 

Figure^ shows how the Su-Olson problem evolves for e = 0.1. Here, one sees how T* (dashed 
lines) and Er (solid lines) approach the same value at t = 250 x 10"^ and the problem reaches 
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thermal equilibrium. Also included in Figure |9] are values for the semi-analytic solution that 
were tabulated by Su & Olson 1996. Circles and squares mark semi-analytic values for Er 
and T"^, respectively. 



Figure [TO] demonstrates how the numerical solution for the Su-Olson problem converges to 
the semi-analytic values. By visually inspecting this plot, one observes two important results. 
First, one notices that the magnitude of the error at a given point decreases by a factor of 2 
as the grid spacing Ax decreases by a factor of 2. This behavior characterizes a first-order 
convergence rate, which one expects for the backward Euler upwinding scheme. Second, this 
convergence test validates that the hybrid Godunov method can capture diffusion behavior 
arising from a system of balance laws with stiff relaxation source terms. When one compares 



the results in Figure 10 to the results in Figure 5 of Reynolds et al 2009, one sees that the 



magnitude of the error is smaller in Figure 10 Calculations were also performed for e = 1 



These solutions also passed through the semi-analytic solutions of Su & Olson 1996. 

10 Full Radiation Hydrodynamics 

One example of full radiation hydrodynamics that has been studied analytically and numer- 
ically is the radiating shock. For sufficiently strong shocks, the post shock material emits 



radiation which penetrates upstream and preheats the system 16,21 . Here, the radiative 
transport of energy plays a significant role in defining the shock structure and evolution [46] . 
Furthermore, ii Pr > p then radiation also transports momentum. 

Lowrie & Rauenzahn 2007 and Lowrie & Edwards 2008 give an excellent summary of the 
solution phenomenology associated with radiating shock waves which are classified into two 
main types: subcritical and supercritical. As the strength of a shock increases, the temper- 
ature behind the shock Ti rises and produces a fiux that penetrates the upstream material 
and heats the region from a reference temperature Tq to temperature Tp > Tq. If Tp < Ti, 
then the shock is termed subcritical. Another characteristic of these shock waves is that at 
higher fiow velocities, Tp increases relative to Ti. At some critical velocity Merit, Tp = Ti and 
such a shock is termed critical. However, if Tp = Ti and u > Merit, then the shock is termed 
supercritical [l6l[2l|[34|[58]. If the material density Pp associated with temperature Tp is 
less than the material density associated with the temperature behind the shock Ti, then the 
material temperature profile has the following features (i) precursor: material is preheated 
ahead of the shock front by radiation, [ii) Zel'dovich spike: overshoot at the shock front 
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that reaches a non-equihbrium temperature T^ > Ti, and {in) relaxation region: as extra 
energy is radiated away, the material temperature declines from Tg to Ti where the width of 



the region is proportional to the post-shock optical depth 34 , 46 



The test problems in this section are based on Ensman 1994, which considered the case of 
a piston moving through static media and followed the evolution with a Lagrangian code. 
Since Nike evolves the dynamical quantities in the Eulerian frame, the problem can be recast 
such that a moving medium impacts a stationary reflecting boundary. The material velocity 
is initialized to match the piston speed specified in the Ensman paper. To produce figures 
that are directly comparable to results obtained from a Lagrangian code, one transforms the 
coordinates of the Eulerian frame calculation according to XEuierian = a^iab — ^inflow^, where 
xiab is the lab frame coordinate, ^Eulerian is the Eulerian frame coordinate, and ^inflow is the 



inflow velocity 21 . The following physical parameters were adapted from Ensman 1994 and 



Hayes & Norman 2003: 



c = 3 X 10"*^" cm s ^, Ooo 



L = 7 X 10^" cm, C = L/Xt 



300 



pO = 7.78 X 10" 



10 g cm-\ T^ 



Physical Parameters: 

3 X 10^ cm s-\ / = 1/3 

22 ^ At = 3.2 X 10^ cm, A^ceii 

T.,oo = {Er^oof^ = 10 K, aa = 3.1 X 10-1° cm-\ 

Subcritical: ^inflow = -6 x 10^ cm s~\ t = [1.7,2.8,3.8] x 10^ s. 

Supercritical: Mmflow = -20 x 10^ cm s"\ t = [4, 7.5, 13] x 10^ s, 

T{x, 0) = Too + 75(x/L) K, Left BC: Reflecting, Right BC: Inflow. 

Motivated by these physical parameters and the non-dimensionalization adopted throughout 



this paper 31 33 34 , one arrives at the following set of values: 



U{x,0) 



m 

E'r 



( 



\ 



P ^inflow 
1„0„,2 I p''T{xfl) 

2r "inflow "•" (7-I) 



r'^(x,o) 



-30/T3(x,0) 



(125) 



Non-dimensionalized Parameters: 

c = lo^ p = 10-", (T„, <Jt = 10, / = 1/3, 
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Figure 11: T and T^ for subcritical shock. 
t = [2,2.5,3] X 10-2. A/'^,u = 512. 
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Figure 12: Fr for subcritical shock, t 
[2,2.5,3] X 10-2. Nc^n = 512. 



At 



uAx 



maxi{\ui\ + acff^ 
T{x, 0) = Too + 7.5 
Subcriticah Mmfiow 



Xin 



0, Xr 



2, AT, 



cell 



512, 



X 



X-n 



-20, t 



-2 



a/4 



[2,2.5,3] X 10" 

Supercriticah MinSow = —66.6, t = [3,6,9] x 10-^, 

Left BC: Reflecting, Right BC: Inflow. 

Figures [TT] and [12] show the results for the subcritical shock. Here, the the Zel'dovich spike 
(temperature overshoot) is sharply peaked, lies just behind the shock, and is only seen in the 
material temperature profile. Additionally, the material and radiation temperatures demon- 
strate non-equilibrium behavior on both sides of the shock front [21]. Furthermore, the flux 
profile is nearly symmetric and matches results discussed by Ensman 1994. As the profile 



moves from left to right, the shock front is located at the point of maximum flux 16 . It is 



important to note how sharply the hybrid Godunov method captures radiative shock fronts, 
where the discontinuity in T is represented with one computational cell. 

Lowrie & Edward 2008 examined semi-analytic solutions of planar radiative shock waves 
with a gray non-equilibrium diffusion model for radiation hydrodynamics and derived an 
estimate of the maximum material temperature at the Zel'dovich spike: 



T. 



LE 



max < Ti 



(3(7.Mg + l)+7P(l-T,^)r 



(126) 
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Figure 13: T and T^ for supercritical 
shock, t = [3, 6, 9] x IQ-^. iVceii = 512. 




Figure 14: Fj. for supercritical shock, t 
[3, 6, 9] X 10-3. iVeell = 512. 



where A^o = ^inflow is the Mach number which determines the strength of the shock. The 
definitions of the other quantities remain unchanged from when they were first presented in 
this paper. Furthermore, the authors cite a coarser estimate of this maximum temperature 
which was originally derived by Mihalas & Mihalas 1984: 



T, 



MM 



(3-7)Ti. 



In Figure [11} the material temperature behind the shock Ti : 
at the height of the Zel'dovich spike T^ ^ 156 for times t 



(127) 

135 and material temperature 
: [2,2.5,3] X 10-2. This value 



compares well with T^ 



LE 
max 



167 as well as the coarser estimate T, 



MM 



180. 



Figures 13 and 14 show the results for the supercritical shock. Here, the shock front is 
located at the temperature maximum. In many instances, computer codes have to resort to 
adaptive mesh refinement to reasonably capture the sharp temperature spike seen in Figure 
13 In these plots, one sees the expected results of an extended region where Tp < Ti and 
the radiation flux Fr being asymmetric |T6}[30|[34]. In Figure 13, Ti = [646,702,738] and 
Tg = [810, 841, 893] for t = [3, 6, 9] x 10-^. For this spatial resolution, these values are within 



the coarser estimate T, 



MM 
max 



[861,936,984]. 



The tests herein examine the behavior of full radiation hydrodynamics. As mentioned ear- 
lier, radiation is the dominant transport mechanism for energy and momentum when Pr > p 
which is equivalent to f'Er > p. If one considers equilibrium diffusion behavior and the 
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Figure 15: A supercritical Zel'dovich spike for various spatial resolutions, t = 0.006. 



material equation of state, then the above inequality becomes fT^ > pT. Since p is initially 
set equal to unity and / = 1/3, this inequality is satisfied by the profiles in Figures 11 and 
T3]- illustrating that radiation can function as the dominant transport mechanism. 



Figure 15 demonstrates self-similar convergent behavior in the numerical solution for a 
Zel'dovich spike in the material temperature T of a supercritical shock wave for a range 
of resolutions A^'ceii = [256,512,1024,2048,4096,8192]. These results were initialized with 
the same conditions as the plot of t = 6 x 10^^ in Figure 13 Figure 15 clearly shows 
asymmetry in the supercritical Zel'dovich spike. This observation corresponds well with the 
discussion in Lowrie & Rauenzahn 2007 and Lowrie & Edwards 2008 as well as the numerical 
results in Sincell et al 1999. However, this physical feature is often missing in many radiation 
hydrodynamical calculations of supercritical shocks. Lastly, the sharp discontinuity from Tg 
to Tp as well as the gradual decay from Tg to Ti in the relaxation region is a testament to 
the hybrid Godunov method in the Nike code. 



11 Conclusions and Future Work 

This paper presents a hybrid Godunov method for full radiation hydrodynamics. Numerical 
tests focusing on the material components, radiation components, and fully coupled system 
show this technique to be uniformly well behaved across various parameter regimes. The 
advantage of this algorithm is that it is accurate, stable, robust, explicit of the material 
flow scale, unsplit, and fully couples matter and radiation without invoking a diffusion-type 
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approximation. Two additional numerical tests that should be investigated are the disper- 
sion analysis of the boundary and initial value problems for radiation hydrodynamical linear 
waves in one spatial dimension. Despite previous computational work being done on the 
boundary value problem by Mihalas & Mihalas 1983 and Turner & Stone 2001, a compu- 
tational test based upon the initial value problem that resembles the propagation of linear 
eigenmodes in hydrodynamics has not been carried out. Constructing such a test seems 
possible given the research of Lowrie et al 1999, Blaes & Socrates 2003, and Johnson 2009. 
These tests are important because they will examine the damping rates and propagation 
speeds of various radiation hydrodynamical waves. 

In this paper, it was shown that the effective eigen-quantities of the modified Godunov 
scheme for the split material subsystem exhibits adiabatic and isothermal behavior in cer- 
tain parameter regimes while evolving the overall radiation hydrodynamical system according 
to an effective CFL condition. Additionally, the modified Godunov scheme resolves the free 
streaming (hydrodynamical) limit with second-order accuracy for the material quantities 
and captures the appropriate material temperature profiles for subcritical and supercritical 
radiating shock waves that are influenced by non-equilibrium diffusion behavior. For these 
reasons, one concludes that the modified Godunov scheme has some asymptotic preserving 
properties. Furthermore, this paper showed that the backward Euler upwinding scheme for 
the split radiation subsystem recovers free streaming, weak equilibrium diffusion, and strong 
equilibrium diffusion behavior even for large cell-optical depths. Additionally, the backward 
Euler upwinding scheme can be cast in a form where its discretization resembles a backward- 
time centered-space differencing of a parabolic/diffusion operator - illustrating consistency 
in the discretization operator. For these reasons, one concludes that the backward Euler 
upwinding scheme has some asymptotic preserving properties. Despite the volume of nu- 
merical tests that are showcased in this paper, it in not clear that the overall algorithm (i.e., 
the hybrid Godunov method which combines the modified Godunov and backward Euler 
upwinding schemes for their respective subsystems) is globally asymptotic preserving. One 
should carry out a discretization analysis on the full hybrid Godunov method in the spirit of 
Lowrie & Morel 2002 which defined in a mathematically rigorous fashion what it means for 
an algorithm to be asymptotic preserving. With respect to this kind of work, full radiation 
hydrodynamics may not be the best system to analytically investigate asymptotic preserv- 
ing properties of the hybrid Godunov method because of the complexity associated with the 
system. One should conduct such an analysis on a set of equations like the one shown in 
Appendix 5 which is a simpler system of hyperbolic balance laws with dual multiscale be- 
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havior. Additionally, it would be interesting to investigate the algorithmic coupling between 
the material and radiation subsystems by implicitly updating the material temperature T 



in a predictor context similar to Equation 124 Such an update could be made after the 



backward Euler upwinding scheme advances the radiation quantities but before the modified 
Godunov scheme advances the material quantities. 

This numerical method employs familiar algorithmic machinery without a significant increase 
in computational overhead. Therefore, it seems reasonable that one will be able to extend 
the ideas in this paper to multiple spatial dimensions via a MUSCL or CTU approach [T]. 
Such work is currently underway. Additionally, the new algorithmic ideas were cast in such a 
way that they should be easily implemented into existing codes, particularly ones that carry 
out MHD calculations. Investigating full radiation MHD in multiple spatial dimensions is 
an area of research on the forefront of computational science and a subject that the authors 
want to pursue. Moreover, one would like to combine a method for radiation hydrodynamics 
with a technique for updating the variable tensor Eddington factor f which is used in the 
closure relation Pr = fE^ for the radiation moment equations. Short characteristic S^ as well 
as Monte Carlo methods are promising techniques for solving the photon transport equation 
and updating f at each temporal iteration (sl. Furthermore, the question remains as to 
whether the hybrid Godunov method and mixed frame approach are applicable to problems 
defined by multigroup physics. Lastly, it would be interesting to compare the performance 
and accuracy of the hybrid Godunov method with an algorithm that is either fully implicit. 



based upon the PnPm family of schemes 11, 12 , or uses a fully relativistic treatment of 
radiation and matter |6 43 . 
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Appendix 1: Modifications to Matrix Equations 
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Outflow Boundary Conditions 
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Reflecting Boundary Conditions 
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Inflow Boundary Conditions 
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Appendix 2: Piecewise Linear Limiting Techniques 

Componentwise Limiting 

1. Compute left, center, and right differences: A_f/j = f/j — f/j„i, Acf/j = ^ (f^i+i — f^j-i), 



2. Enforce the TVD condition: An^M = 2min (|A_f/i|, \A+Ui\) 
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3. Define Hmited slope: 

P (AU) = l ™i'^(l^-^^l' |Aiimf/i|)sign(A,t/i) A_t/,A+f/, > 
'^ \ A_UiA+U^ < 

Limiting Across Characteristic Fields 

1. Compute left, center, and right differences like Step 1 above 



2. Using the left material eigenvectors L^ given in Equation 78, expand the differences 
in characteristic variables where k defines a specific row vector of L^: AJJ^ = L^' ■ 
A_f/„ A,Wf = L^/ ■ A,f/„ A+Wf = L^/ ■ A+f/, 

3. Enforce the TVD condition one field at a time: Aii^Wf = 2min (|A_Wf |, |A+Wf |) 

4. Define limited characteristic variables: 

min ( I AMt I , I AiimWf | ) sign ( A^Wf ) AMfA+Uf > 



AU 



A_WAM^ < 



5. Using the right material eigenvectors R^ given in Equation 77, reconstruct the limited 
slopes for each of the material quantities where k defines a specific column vector of 

Appendix 3: Material Flux Evaluations 

The Lax-Priedrichs flux is first-order accurate for smooth solutions, where: 



1 Ax 



FZi2 = 7T {FiUr.) + F{Un)) + -— (f/^ - Ur) . (128) 



The Richtmyer fiux is second-order accurate for smooth solutions and is based on an 
intermediate quantity Uj]^^,^, where: 

^^1/2 = nulii2). Uli/2 = \^^^ + ^r) + \^^ (^(^^) - ^(f^«)) • (128) 

The HLLE fiux is higher-order accurate for smooth solutions and is derived by approximately 
solving for the Rankine-Hugoniot shock conditions of the Riemann problem, where: 



F{Ul) < SL 

1/2 - 



i^HLLE ^ ^HLLE ,^ < Q < .^ , (128) 



F{Ur) Sr < 
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nHLLE 



srF{Ul) - slF{Ur) + slSr {Un - Ul) 



(128) 
sr- Sl 

Here, sl,r are the wave speeds to the left and right of a cell interface that can be estimated 

in the following three ways in the algorithm (?) directly: Sl,r = Ul,r =F aeflf(f^L,_R), (ii) using 

Roe-averaged quantities: sl^r = ul,r =F acff(f/L,i?) where R = {p/ kl + Pr k.r)/{Pl + Pr ), 

and {in) taking the maximum wave speed over all grid cells: s^^r = =FSmax where Smax = 

maXj(|-Uj| + acs,i)- It is important to note that when sl^^i = =FSmax the HLLE flux function 

algebraically becomes the Lax-Friedrichs flux function. 



Appendix 4: Boundary Conditions for Corrector Scheme 

If one assumes that there is a computational grid of A^ cells starting at the index istart and 
ending at the index i^nd and if one further assumes that there are Ng^ost cells on either side 
of A^, then one defines boundary conditions in the following ways: 



Periodic Boundary Conditions 



Left: 



Right: U{i 



Uyistart — Nghosti • • • J istart ^ 1) 



end 



■ 1 ^end 



N. 



ghost 



U{iend — Nghost + 1, • • • , lend), 
^ \^start7 • • • 7 ^start ' ^^ ghost -*-/• 



Outflow Boundary Conditions 



Left: U{istart — Nghost, ■ ■ ■ , istart — 1) ^ U{istart), 

Right: U{iend + 1, • • • , iend + Nghost) -^ U(iend)- 



Reflecting Boundary Conditions 

Left: {p, E, Er}{istart — Nghost, ■ ■ ■ , istart — 1) 

\Tn, ^riy^start ^^ ghost, • ■ • , ^ start ^) 

Right: {p, E, Er}{iend + 1, • • • , iend + Nghost) 
{m, Fr}{iend + 1, • ■ • , iend + Nghost) 



{P, E, Er}{istart + Nghost ~ 1, • • • , istart), 
\VTi, i^^rix^start ' ^^ ghost ^, • • • , ^start), 
{P, E, Er}{iend, ■ ■ ■ , iend — Nghost + 1), 
— {m,Fr}{iend, ■ ■ ■ ,iend — Nghost + !)■ 



Inflow Boundary Conditions 

Left: Ud.fnrt — N, 



^ \^ start ^^ ghost , • ■ • , ^ start -*- 1 
Right: U{iend + 1, • ■ • , iend + Nghost) 



^ x'^start ^^ ghost-; ■ • ■ i '^start ^J-) 
U^iiiejid + 1, . . . , iend + Nghost)- 



A Hybrid Godunov Method for Radiation Hydrodynamics 



61 



Appendix 5: Simpler System for Investigation 

To illustrate how some systems of hyperbolic balance laws with dual multiscale behavior 
reduce to non-stiff systems, consider the following set of partial differential equations: 

(125) 
(126) 
(127) 
(128) 

where the terms appearing in the above system do not directly correspond with material and 



Pt + m^ = 


= 0, 


mt + p^ = 


= ^/, 


et + Cf^ -- 


= Ca{p-e 


ft + Ce, -- 


= -Caf, 



radiation quantities. Equations 125 and 126 define the slow (i.e., macroscopic) hyperbolic 



system, which has a wave speed of unity; while Equations 127| and 128 define the fast (i.e., 
microscopic) hyperbolic system, which has a wave speed of C. Stiffness associated with 
the different wave speeds defines one multiscale behavior. Additionally, the non-zero source 



terms in Equations |126p28| are augmented by a magnitude a which may be small or large 
with respect to unity and defines the other multiscale behavior. In some parameter regime 
of C and a, one expects e — )■ p because the source terms enable the relaxation of some of the 
components in the system. If one assumes that there is a small parameter e ^ 1 such that 
C,a — !■ 1/e and one employs methods from perturbation theory to construct the following 
series for e and /: 



e = Co + eei -F e^e2 + e^Cs + . . . , 



f 



/o + e/i + eV2 + e'/s + 



(129) 
(130) 



after matching terms of like order in e, one arrives at the following reduced system of equa- 
tions which exhibits no stiffness and has a well-defined solution: 



Pt + m^ = 0, 
mt + 2p^ = + C(e2), 



(131) 
(132) 
(133) 



